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We numerically calculate the dissipative part of the self-force on a scalar charge moving on a 
circular, geodesic, equatorial orbit in Kerr spacetime. The solution to the scalar held equation is 
computed by separating variables and is expressed as a mode sum over radial and angular modes. 
The force is then computed in two ways: a direct, instantaneous force calculation which uses the 
half-retarded-minus-half-advanced held, and an indirect method which uses the energy and angular 
momentum flux at the horizon and at infinity to infer the force. We are able to show numerically 

■ and analytically that the force-per-mode is the same for both methods. To enforce the boundary 

conditions (ingoing radiation at the horizon and outgoing radiation at infinity for the retarded 
solution) numerical solutions to the radial equation are matched to asymptotic expansions for the 
fields at the boundaries. Recursion relations for the coefficients in the asymptotic expansions are 

' given in an appendix. 
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INTRODUCTION: MOTIVATION AND SUMMARY 



With sensitive ground-based gravitational wave detectors, such as LIGO, GEO, TAMA and VIRGO 0, H m 
operation and the launch of a space-based gravitational wave detector LISA planned in the next decade, the need 
to accurately model the waves emitted from inspiraling binary star systems has become acute. 

In the case of LIGO and the other ground-based detectors, the sensitivity of the instrument is at a maximum around 
a few hundred Hertz. This frequency is comparable to the orbital frequency of two neutron stars just prior to their 
, final coalescence, and therefore LIGO is particularly sensitive to these sources. These constituents of these systems 
have comparable mass and relatively weak gravitational fields; therefore the emitted signal can be computed using a 
\ weak-field, slow-motion approximation (i.e a post-Newtonian approximation Although this method - and other 
improvements to the post-Newtonian method |fj - give waveforms that are accurate enough to detect signals buried 
| in the LIGO noise, a more detailed knowledge of late inspiral and merger, accessible only by numerical evolution, is 
needed to extract from the waves the astrophysics of their sources. 

LISA, however, will be sensitive to signals with much lower frequency, in particular, to waves from a stellar-size black 
(3JT), hole spiraling in to a (super-)massive black hole [10 3 — 1O 6 M ]. Detectable sources with this extreme mass-ratio will 
have periods of (many) minutes and will persist for weeks (or even months). In addition, the smaller mass may spend 
part of its orbit deep in the strong gravitational field of the larger mass, making a post-Newtonian approximation 
inappropriate for computing the predicted waveforms. Clearly, this problem is better suited to black hole perturbation 
theory where the background geometry generated by the larger mass M is treated exactly, and the smaller mass fi 
generates a small perturbation of the geometry. To lowest nontrivial order in fi/M, the field generated by the smaller 
mass causes its trajectory to deviate from a geodesic of the background spacetime. This is the origin of the self-force. 
Hcuristically, one can think of the smaller mass as traveling on a geodesic of the perturbed spacetime. Although 
this description of the force is intuitive and compelling, computing the self-force entails a number of conceptual and 
technical difficulties. 

At the center of the conceptual difficulties is the renormalization problem: even though the perturbing mass is 
small compared to the central mass, the perturbed field diverges at the position of the particle. A number of authors 
0, IS 0> EH EH E3 have addressed this difficulty and developed formal methods for removing the divergence while 
preserving the finite parts of the field that give rise to the self-force. These papers solve the problem in principle; there 
remain, however, a number of difficulties in implementing these prescriptions to find the force and subsequent motion 
and waveform. In particular, there remains a gauge problem. The self-force equations for a small mass moving in a 
background spacetime have only been written down in Lorentz gauge, but, unfortunately, most methods for finding 
the perturbing gravitational fields use a different gauge. Complicating matters further, recent calculations [l3l | suggest 
that the gauge transformations that relate gauges of perturbation theory to Lorentz gauge are poorly behaved. Thus 
it is difficult to use the metric perturbations that are readily available from black hole perturbation theory in the 
formal equations for the self-force. 

One method of side-stepping both the gauge and the renormalization issue is to compute the self-force indirectly 
by computing the energy and angular momentum flux at infinity. Such techniques usually assume that the particle 
is in circular geodesic orbit. The particle is then the source of the metric perturbation (or of the electromagnetic or 
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scalar field perturbation). The perturbed field is examined at the horizon and infinity and the energy and angular 
momentum flux is computed. The rate of energy loss is then equated to an orbital energy loss and the rate of inspiral 
can be inferred. Although such techniques have been broadly applied (e.g. 013) they have a drawback: They give 
only the time-averaged dissipative force. [See ^(| for a discussion of the shortcomings of the time-averaged force.] 
When used to infer the rate of orbital decay, these energy balance arguments assume that the field energy on a spatial 
slice is constant as we go from one orbit to the other. See Section IVA for discussion. 

In this paper, we attack a somewhat simpler problem that still puts us on a path to solve the more general problem. 
We assume that our particle is a scalar charge that acts as the source of a Klein-Gordon field. The self-force arising 
from the back-reaction of a scalar (rather than a gravitational) field involves no delicate gauge issues. The part of 
our work reported here avoids the renormalization problem by using the half-retarded-minus-half-advanced field. The 
divergent structure of the field at the position of the particle is the same for the advanced and retarded solutions, and 
therefore the difference is smooth at the particle's position. By evaluating the gradient of the half-retarded-minus- 
half-advanced field we are able to compute the dissipative (time-antisymmetric) parts of the self-force acting on the 
particle at any instant; our method does not require time-averaging to compute the force. Although in the present 
paper we will assume that our source particle is traveling on a circular geodesic orbit, it is not an essential feature 
of the method. Extending the method to a particle in non-circular orbit would severely complicate the numerical 
calculation of the field, but it would not appreciably complicate the calculation of the self-force. However, in this 
special case of circular motion, we are able to show that our direct calculation of the instantaneous force agrees - 
mode-by-mode - with the indirect calculation based on the average energy loss at infinity and at the horizon. 

The next steps in this work are to perform calculations for non-circular orbits and to calculate the conservative as 
well as the dissipative part of the self-force. As discussed above, the only complication arising in the former task is the 
numerical calculation of the field from an orbit that has a countably infinite set of frequencies. The latter, however, 
requires a new procedure that will include a renormalization of the divergent field. It is in anticipation of the delicate 
subtraction arising in mode-by-mode renormalization that we have taken such care with the accuracy of the radial 
functions computed in this paper. 

Outline and conventions 

In section [H] we use separation of variables to solve the wave equation and write the scalar field $ as a sum over 
radial and angular mode functions. In section ITTT1 we describe the scheme for the numerical construction of the radial 
functions which appear in this solution. In section HVI we present the argument that the two methods are equivalent, 
and derive separate expressions for the evolution of the conserved quantities E and L based on each one. Finally 
in section we discuss numerical details, show the properties of the radial functions, and present our results for 
dissipative self-force. 

Throughout the paper we use Boyer-Lindquist coordinates (t,r,6,(f>) and a metric with signature ( — h H — h). The 
mass of the black hole is M, its spin is a, and the usual abbreviations A = r 2 + a 2 — 2Mr and £ = r 2 + a 2 cos 6 
are used; the radial coordinate of the horizon is written r + = M + V M 2 — a 2 . Our conventions are those of Misner, 
Thorne and Wheeler ^3] > except that their symbol for £ is p 2 . 



II. SOLUTION OF THE WAVE EQUATION 



We consider the Klein-Gordon field $ of a point paticle of conserved scalar charge q, orbiting on a circular, equatorial 
geodesic of the Kerr geometry. The field $ satisfies the massless scalar wave equation 

V Q V Q $ = -4ttp, (2.1) 

with the scalar charge density p a delta function along the trajectory z°(t). 

P = qJs 4 (x a ,z a (T))dr. (2.2) 

With this normalization, we have q = / pdV, where dV is the volume element on a spacelike hypersurface V 

Jv 

orthogonal to the trajectory. For a circular orbit of radius ro and frequency Q, we have 



P = Q ^-t^-t — dr, (2.3) 

J £ sin 

where it* is the (constant) time component of the particle four- velocity, and Ea. l|2.I|l takes the form, 

V a V Q $ = -An-?— Sir - r )8(9 - VjSU - tit). (2.4) 
riu* 
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Separating variables yields the solution 

$ = ^■ S tY,Ri m (r)S tm (z)S lm (e)e im «- at \ (2.5) 

1,771 

where the Si m and Ri m satisfy angular and radial equations given below. The Si m are oblate spheroidal harmonics, 
given by 



1 



sin 9 89 



d ( 8 
sm 9— 



sin 2 9 



m 2 n 2 a 2 cos 2 9 



Sim (9) = 0. 



We take the harmonics to be real, and normalize by 

S lm (9) 2 sm(9)d9 = 1. 



(2.6) 



(2.7) 



We fix the sign by demanding that Si m e %m ^ = \/2irYi m when a = (A; TO will reduce to 1(1 + 1) in this case). The 
numerical calculation of the Si m and their eigenvalues A; m is a straightforward task, and is not discussed here. (See 
Numerical Recipes [T^| for a treatment similar to ours.) The computation of the Ri m is considerably more complicated 
and is detailed below. 



III. INTEGRATION OF THE RADIAL EQUATION 

The radial functions satisfy 

8 ( . d \ m 2 /„, , 9 9N 2 



^Uj-j+^fn 2 (r 2 +a 2 ) 2 - iQMar + a 2 ) - m 2 Q 2 a 2 - X lm 



Ritu (r) = -iwS (r - ro) ■ (3.1) 



First found by Carter [T<|, this is also the Teukolsky equation with spin s = 0, specialized to our source (this equation 
can be derived from Teukolsky ( 20] ) by setting u) — mVL). To provide accurate initial data for numerical integration 
of this equation, we solve the equation in asymptotic series valid near the horizon and near infinity. In terms of the 

angular velocity of the horizon, uj + = , and the "tortoise coordinate" r*, satisfying 

2Mr+ 

the series take the form 

Kn = E °n (r r+) n e - im ^-»+>* (near r + ), (3.3) 

R ^ = E it emnr * ( near °°)- ( 3 - 4 ) 

n=l 

The signs in the exponentials of these series amount to a choice of boundary conditions. Here we have picked 
ingoing radiation (— ) at the horizon and outgoing radiation (+) at infinity, to construct the retarded solution. Note 
that in the case of circular orbits, where the time dependence is just e _ ™ a *, observers at infinity and the horizon 
always agree on the direction of the radiation down the hole, and (— ) is always the correct sign for this solution. 

The coefficients and satisfy recursion relations given in the appendix, with overall normalization set for 
each I, m harmonic by the <5-function source. The asymptotic series provide initial values of Ri m and its derivative 
for numerical integration of the homogeneous radial equation out from the horizon and in from infinity. We thereby 
obtain homogeneous solutions R^ (from the horizon to ro) and Rf^ (from ro to infinity). To obtain the solution Ri m 
to the inhomogeneous wave equation (|3.1|) . we patch R^ m and RfZ at r = r$. Requiring that Ri m be continuous and 
that the discontinuity in dRi m /dr be fixed by the (^-function source, we have 

Rl m = C lrn R+ m (r < )R% l (r > ), (3.5) 
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where r< (r>) denotes the lesser (greater) of r and r, and 

— 47T 

r- (3-6) 

A-W(R+ m ,RZ) 

It is easy to show that the quantity in the denominator, A times the Wronskian of the integrated solutions, is a 
constant for any two solutions of the homogeneous radial equation. It will be convenient to define constant factors 
Cf m and Cj^ for which R im = C^ m R^ m when r < r and R lm = Cf^Rf^ for r > r , namely 

r + = -47rflj%(r ) 
lm AW(R+ m ,RZ) { ' ' 

Clm Aw(R+ m ,Rr m y 



IV. RADIATION REACTION 
A. Equivalence of the Methods 

The self- force on a body can be split into two pieces, a conservative force and a dissipative force, by writing the 
retarded field of the source as a sum of parts that are even and odd under the interchange, advanced <-» retarded, of 
ingoing and outgoing radiation. 

$rct = ^($rot + $adv) + ^($rct - $adv) (4.1) 

(Because advanced and retarded fields are well-defined on generic, globally hyperbolic spacetimes, this division is 
generic for the self- force arising from linear fields.) The self- force that we compute is linear in the field $ and can 
similarly be written as a sum of a part invariant under the interchange advanced <-> retarded - the conservative part 
of the force, and a part that changes sign under the interchange - the dissipative part of the force. 

For a charge sustained in perpetual circular orbit by an external force, the field energy is unchanged from one 
hypersurface to the next, and the work done by the self-force is therefore equal to the energy radiated to the horizon 
and to infinity. Because the sign of the radiated energy changes under the advanced <-» retarded interchange, so does 
the sign of the work done by the self-force. Thus the work done by the self-force is the work done by its dissipative 
part. 

Formally, the stress tensor of a finite test-mass fj, is the sum of contributions from the mass, the scalar field, and 
the external force, 

rpaj3 rpaf) _|_ jia/3 . rpOtff 

with 

T f = ^ - i 9ct/3 V 7 $V^$) . (4.2) 

Associated with the timclikc Killing vector t a is the conserved current j a — ~T a ^ . The self- force is constructed 
from —\7pTg 13 = pV"^, and the work done by the self-force between two t — constant hypersurfaces Si and £2 is 

W= f d 4 Vpt a V a $ = f d 4 VV aJ %, (4.3) 
Jn Jn 

with f2 the region of spacetime between £1 and £2- We can write this integral as the sum of a vanishing time derivative 
and the integral of a 3-dimensional divergence by writing 

V Q . ?s tt = d t f + ^D a {Nf), (4.4) 

where TV is the lapse, D a the covariant derivative operator on the hypersurface, and jg the projection of jg into the 
hypersurface. Using 

/ ddWV =( f - f ) .f s Nd 3 V - 0, 



5 



we have 

W = / j%dS a + [ j%NdS a , (4.5) 

JSoo iff 

where, by J s is meant the limit linv^oo J„ , with SV a sphere of constant coordinate r; and -ff is the horizon. 
Because t a is orthogonal to the horizon and to S r , W has the form 

An 



D a <f>d t <S>dS a + / L> Q $<9 t $iVdS a 



The change of sign of these flux integrals under the advanced <-> retarded interchange, clear on general grounds, is 
apparent from this last form, together with the fact that the interchange is equivalent in the Kerr geometry to the 
transformation induced by the diffeo t,<f> «-> — t, —cj). The explicit form of these integrals for each mode is given in the 
next section. 

The self-force measured is balanced here by an external force. When no external force is present, and when one can 
neglect the change in the radiative part of the field energy on successive hypersurfaces, the work W is equal to the 
change in the energy of the particle between the hypersurfaces Si and £2. Modeling the particle by a family of dust 
balls with stress-energy pu a u@, with it Q u*(i Q + Q4> a ), we have 



W = / V a T^t l3 d i V= / V a ( P u a u f3 t f3 )d 4 V 

u a V a (pu^)d A V. (4.6) 
If, along a circular geodesic in the mass, we define a momentum for which 

Pa t a = I TfotPdSa = I puptP^NdV, (4.7) 

JT, 

/d 
— (pBt /3 )Ndt, and we have 
dr 

Because the dissipative field, S'diss := |(^ret — ^adv), is regular for a point particle, the dissipative part of the 
self-force is well-defined for a point particle, without renormalization. The work done by the self-force at the particle, 
must then, in our case of perpetual circular motion, be identical to the radiative flux of Ea. (|4.5|l . 



B. Energy and Angular Momentum Flux 

We have just argued that we can indirectly find the self-force on the charge by computing the flux integral of 
Eq. (|4.5|) . The flux of angular momentum to infinity and to the black hole is similarly given by the integral 

/ IsdSa + I mdS a , (4.9) 

with jg = TpcjJ 3 , with <jf the rotational Killing vector 9^. 

Again, when no external torque is present and when the difference between the field angular momentum on successive 
hypersurfaces can be neglected, the radiated angular momentum is equal to the change in the angular momentum of 
the particle. 

We denote by E — —u t = u a t a and L = = u a (j) a the energy and (z component of) angular momentum per unit 
rest mass of the charge. In the full dynamical problem, no flux reaches spatial infinity, and one relates a change in 
energy between successive hypersurfaces to asymptotic flux by choosing a family of asymptotically null hypersurfaces. 
In our model, however, the charge has been orbiting forever, and the integrals can be evaluated at spatial infinity and 
at the bifurcation horizon (as in the previous section). Consider first the flux of angular momentum L. We compute 
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Ea. H4.9|l to get (times the rest mass /z), for the particle. For a surface of constant r (and t), the flux in the direction 
of increasing r is 

C = J T^EtKl = J T r4 ,Adn. (4.10) 

The real scalar field is a sum $ = ^ $; m of complex fields whose angular eigenfunctions Si m e tm ^ are orthogonal. 

The flux integral is then sum of integrals for each mode, and each mode integral involves a stress-energy tensor of the 
form 

T *P = ^ {^{A^^lm - 5^V 7 #£ B V 7 ^Im) • (4-11) 

The modes have asymptotic behavior 

$im = 2^°^ mSlm ® Slm {6) e- im( - n -"+>* j^-W + 0(r - r+) (near r+) (4.12) 

= AC^m (f) (^) -e lmS!r *e tm W- m ' + 0(r- 2 ) (near oo). (4.13) 

Ziru 1 r 

The (r — (ft) component of the stress energy in each mode has corresponding behavior 



( T rAm = -^- r -^-nm 2 <£ lm [1 + Oir- 1 )] (near oo), (4.15) 



47T A 

and by Ij4.1t)jl the fluxes are 

£M = _ ^(2^) 2 " m2|C ^ |2 ^ m( ^ 2 - (4 - 1?) 

These give the rate at which angular momentum crosses the surfaces in the direction of increasing r. The flux 

integrals (j4.5jl , 14.9jl involve the outward normal of a 3- volume bounded by spatial infinity and the horizon, implying 
a sign change for the contribution (|4.16|) at the horizon: 

~' r $> 2 S im (f ) 2 \2Mr + (fi - u + ) \C+J + \C^f\ . (4.18) 



dt 87T 2 /i(^* 

For 2£, a similar calculation yields 



u i , 



dE dL 



This relationship is in fact necessary for the particle to remain in circular orbit as it radiates and spirals inward. 

C. Half- Retarded Minus Half- Advanced 

Here we calculate the dissipative self-force directly, by taking the gradient of the half-retarded-minus-half-advanced 
solution. The expression of (j2.5jl is the retarded solution. As noted above, the advanced solution is obtained from the 
retarded by the transformation (ft —> —(ft, t — ► — t. Denoting the retarded solution by (— ) and the advanced solution 
by (+), we have 

*± = iE^^^^ 1 *" " 1 (4 - 20) 

l,m 
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and the dissipative solution is then 

$dis = \ - $+) = ^j^ii2i m (r)S, m ('/2)Sim(e)Bin[m(^- fit)] . (4.21) 

To elucidate this expression, note the m — > —m symmetry properties of its components. Eqs. (|4.12() and 14.13|) show 
that Ru- m ) — R* m - The Si m were chosen to have the same sign as the Legendre functions P™, so Sj(_ m ) = (— l) m Si m . 
The sine function of course changes sign with m. This all implies that the imaginary part of each mode of <I>dis (coming 
from the real part of Ri m ) changes sign with m while the real part doesn't. Since m = does not contribute to the 
sum by virtue of the sine, the solution is real (like the retarded and advanced solutions), and more importantly comes 
completely from the imaginary parts of the radial functions, which, as solutions to the homogeneous radial equation, 
are smooth. We thus obtain the more transparent version of <&disi 

$ dis = Yl Im [^™( r )] Si m (*/2)Si m (9) sin [m (0 - fit)] , (4.22) 

Lm 

each term of which is smooth at the particle. The force is given by minus the orthogonal projection of the gradient 
of this quantity (evaluated at the location of the particle), times the scalar charge q. To wit, 

(p a ) = -q(5i + U a U )V p Q Aia . (4.23) 

ar 

The location of the particle is (j> = fit, so the r and 9 components vanish by virtue of the sine in Eq. I|4.21[l . The t 
and cj> components give us the evolution of the conserved quantities E and L. Because V* = — fW^, the projection 
operator is just the identity (ie, (6% + u Q u /3 )V/3 = V a ), and Ea. i|4.23|) yields 

AT _„2 

E mlm i R i™ M]S* m (V2) 2 (4.24) 



dt 2tt^ (w*) z 



Lm 



We have used 4- — u*^ to rewrite for easy comparison with Ea. (|4.18|l . Eqs. (|4.18l) and l|4.24|l must be equal, and 
furthermore the equality must be mode by mode, since the argument of section HV Al is valid mode by mode. We can 
therefore write 

47rIm[fl, m (ro)] =m[2Mr + (n-Lu+)\C+ m | 2 +ft|C^| 2 ] , (4.26) 

and in fact this relationship provides a convenient means to estimate the total numerical error involved in the com- 
putation of a given Ri m . 

V. DISCUSSION 
A. The Radial Functions 

1. Computation 

The general scheme for the computation of the radial functions is presented in Section II I II To summarize, one 
integrates in from infinity and out from the horizon, using asymptotic series to provide initial values, and then 
patches the two solutions together at the particle. Here we discuss our actual implementation of this method. There 
are two sources of error in the computation: the integration itself, and the initial value that starts it (gotten from the 
series). We pick a desired maximum error for the whole computation and make sure that the error from neither source 
exceeds this bound. For the integration, we use the Runga-Kutta Cash-Karp order 4/5 routine, as implemented in 
the GSL scientific library |2l| . which reports the error introduced at each step (note that we integrate in r, not r*). 
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The sum of the magnitudes of these errors is a (gross) upper bound on the total numerical error, and we adjust the 
parameters of the integrator until this falls below our desired maximum. We control the error introduced by the 
initial value by increasing the value's accuracy until the result of the integration (at ro) does not change to within 
the desired error. One can increase the accuracy either by increasing the number of terms in the series used or by 
moving further away from the particle, where the series are increasingly convergent. We chose to fix the number of 
terms and increase the distance. For the results accurate to six significant figures presented in this paper, we used 
fifteen terms from each series, allowing us to begin most infinity integrations somewhere between r/M = 100 and 
r/M = 1000, and most horizon integrations at distances of order 10 -2 to 10 -3 away from the horizon. The series 
appear to converge slowest (i.e., we must begin furthest from the particle) when the difference between I and \m\ is 
largest. The patching of the solutions requires the calculation of the Wronskian; although this can be done at the 
single point r , we integrate a little past r in either direction to have a range of values over which we can verify that 
the quantity AW(R/ nl , R^) of Eci. (|3.6|) is constant (to the desired accuracy). Finally, at the end of each integration 
we check that Ea. (|4.26|l is satisfied (again to the desired accuracy). 



2. Properties 

Although the field they sum to is real, the radial functions themselves are complex-valued. At the particle, the 
real parts have cusps while the imaginary parts remain smooth (see Fig. ^for a picture). This is because the delta 
function in our source is purely real. The imaginary parts can be thought of as solutions to the source-free wave 
equation, coupled to our problem via the radiative boundary conditions. Only these smooth functions contribute to 
the dissipative field (see Eq. I4.22[l . Away from the particle, the radial functions can be understood in terms of the 
asymptotic series Eqs. (|3.3fl an d (|3.4|l . The first terms of these series describe the end behavior, 



Rim (r -» r+) ~ e - im ^+>*, 

R lm (r^cx>)~ i e imQr '. 

r 



(5.1) 
(5.2) 



To understand these we must 'translate' from r* to r. As r — ► oo, r* — ► r, so the behavior at infinity is simply a 1/r 
decaying sinusoid of frequency mQ. As r — > r + , however, r* — > — oo. This means that the distance between points in 
r is increasingly stretched by r* as one approaches the horizon. The distance between wave crests is constant in r*, 
so in r it is steadily decreasing. The effect is that wavefronts 'pile up' on the horizon, as in Fig. 2(a) below. 
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FIG. 1: Two Radial Functions, for a = .8M, ro = 7M (prograde orbit, w .05). The real parts have cusps at the position of 
the particle, while the imaginary parts are smooth everywhere. The oscillations in (b) are too small compared to the cusp to 
be seen (we show these below). 
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(a) Near the horizon, = 1.6M 



(b) Seeing the infinity behavior 



FIG. 2: Closer looks at -Rio.io of fig. |l(b)| At the horizon wave fronts pile up, and at infinity we have a decaying sinusoid. 



B. Dissipative Self-force 

The dissipative self- force is characterized by calculated from Eg. 1)4. 18(1 or equivalently Ea. (|4.24[) . A particle 
in the adiabatic regime will inspiral along circular orbits according to this quantity until it begins a transition to a 
'plunge' orbit near the inner- most stable circular orbit (ISCO), as discussed in [23. Thus it is only meaningful to 
compute results for orbits outside the ISCO. The convergence of the l,m sums that give ^ is quite good, until the 
particle gets too close to the black hole. For most cases, the ISCO is far enough out that all orbits of interest require 
a small number (I = 1Q — 15) of modes to converge to six significant figures. However, as a co-rotating hole moves 
the ISCO inwards (see |23j for a nice discussion of circular orbits in Kerr), more modes can be required. For the 
extremal co-rotating case, the ISCO approaches the horizon at r = M and the adiabatic regime involves orbits whose 
convergence is quite bad. For the r = 2M orbit result shown in fig. |31 below, I = 30 modes were required for the six 
significant figures. 
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FIG. 3: Some values of % x 4 for various orbits. Negative values of a represent retrograde orbits. Results for orbits within the 
inner-most stable circular orbit are unphysical and not shown. All results are accurate to the six significant figures displayed. 
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VI. APPENDIX— THE RECURSION RELATIONS 



The C+ and of Eqs. (|3.3|1 and (|3.4I) satisfy recursion relations, which were derived with the help of Mathematica. 
The procedure is to plug in the forms of Eqs. (|3.3|l and l|3.4|l to the radial equation (|3.1|l . and rewrite the results with 
a single summation index, so that one derives an expression that must be zero (the equation is homogeneous away 
from the particle). These expressions are 

(6.1) 
0, (6.2) 
with ft and f°° given below (where we have defined Q+ = O — oj+). 



i=0 
5 

ErCX) s~lOO 
J i ^n—i 



f+ = -Aam 2 Mr+n + 4:M 2 [n 2 -2imnr + n + + 2m 2 Mr + (n 2 -n 2 + )] 

+a 2 [-An 2 + 8imMnfl+ + m 2 (l - 4M 2 (fl 2 - Q 2 + ))] 
/+ = -2r+ (-1 + 3n — 2n 2 + A + a 2 m 2 fl 2 ) + Aia 2 m (-1 + 2n) Q+ - 16m 2 M 2 r+ (fl 2 - fl+ 2 ) 

+2M [l + 2n 2 - A + 2am 2 n + 3a 2 m 2 D, 2 - Aimr + fl + - Aa 2 m 2 Q + 2 + n (-3 + 6imr+fi+)] 
/+ = 2 + n 2 - A - ha 2 m 2 tt 2 + \2m 2 Mr+n 2 - AimMQ. + + 10imr + n + + Aa 2 m 2 n + 2 

-12m 2 Mr + ^+ 2 + in [3i + 4m (M - 2r + ) Q + ] 
f+ = 2m [-i (-2 + n) fl+ + 2mr + (fl 2 - fl+ 2 )] 

ft = m 2 {n 2 - 0+ 2 ) 

/ °° = -2im(-l + n)Q 

f^° = 2 + n 2 - A - 8imMfl - a 2 m 2 fl 2 + n (-3 + UmMfl) 

/ 2 °° = 2[ia 2 m(5-2n)n + M (-10 + 9n-2n 2 + \-2am 2 n + a 2 m 2 n 2 )] 

/ 3 °° = 4M 2 (-3 + nf - a 4 m 2 n 2 + a 2 [18 + m 2 - I2n + 2n 2 — A + AimM (-3 + n) Q] 
/ 4 °° = 2a 2 (-4 + n) [M (7 - 2n) - ia 2 mVL\ 
/ 5 °° = a 4 (-5 + n)(-4 + n) 
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